(N 

o 

(N 
X 



Measures and LMI for space launcher 
robust control validation* 

Didier Henrion 1 ' 2 ' 3 , Martine Ganet-Schoeller 4 , Samir Bennani 5 

March 3, 2013 



Abstract 

We describe a new temporal verification framework for safety and robustness 
analysis of nonlinear control laws, our target application being a space launcher 
vehicle. Robustness analysis, formulated as a nonconvex nonlinear optimization 
r~| ' problem on admissible trajectories corresponding to piecewise polynomial dynam- 

ics, is relaxed into a convex linear programming problem on measures. This infinite- 
dimensional problem is then formulated as a generalized moment problem, which 
allows for a numerical solution via a hierarchy of linear matrix inequality relax- 
ations solved by semidefinite programming. The approach is illustrated on space 
launcher vehicle benchmark problems, in the presence of closed-loop nonlinearities 
(saturations and dead-zones) and axis coupling. 
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1 Introduction 

X 

This work is carried out within the scope of the project SAFE-V (Space Application 
Flight control Enhancement of Validation Framework). The objective of this project is to 
analyse, develop and demonstrate effective design, verification and validation strategies 
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and metrics for advanced guidance navigation and control systems. These new strategies 
should be implemented in an incremental manner in the traditional validation framework 
in order to be applicable both for current launchers validation and for future launcher 
and re-entry vehicle validation. 

Control design performance and compliance with system specifications should be demon- 
strated by a suitable combination of stability analysis and simulation. This validation is 
made on a wide domain of variation of influent parameters; it includ es worst case vali- 
dation and Monte Carlo simulation for statistic requirements, see e.g. iRongier and Droz 
(119991 ) for Ariane 5 validation. 

Generally speaking, the Monte Carlo method is a numerical integration method using sam- 
pling, which can be used, for example, to determine the quantities of interest of a variable 
of interest for a random input variable such as mean or standard deviation or probability 
density function. One of the most commonly used methods inside ASTRIUM ST for th e 
confidence interval on a quantile estimation is Wilks formula, see lOPENTURNSl ((2008). 
This formula of the confidence interval depending on the probability level and the number 
of samples can be used, for Monte Carlo only, in two ways: to determine the probability 
and confidence level for the values of samples chosen by the user, or in reverse to determine 
the number of simulations to be carried out for the values probability and confidence level 
chosen by the user. Monte Carlo simulations are widely used in traditional validation 
process and they are very robust because they do not require smoothing assumptions on 
the model neither importance factor assumptions on the parameters. The random events 
are simulated as they would occur naturally. The generic Monte Carlo method is however 
very time consuming, because a great number of simulations to demonstrate extreme level 
of probability (in general much more than 1/p simulations to estimate the probability of 
an event of probability p). The number of simulations required may be critical for long 
duration mis s ion si mulation (like launcher flight). In the context of robust control, see 
Tempo et al.l ( 120051 ) for related probabilistic approaches. Note that these approaches may 
miss worst case behavior, especially when the number of parameters is large. 

In contrast, worst case validation techniques are aimed at identifying worst case instances 
and behavior. These techniques include Lyapunov and mu-analysis approaches. Among 
them, the most relevant for space vehicle applications are: 



results by the GARTEUR g roup, see Fielding et al.l ( 120021 ). and the COFCLUO 
group, see ICOFCLUOl (120071 ). for recent aeronautic applications and new develop- 
ments; 

results by AST RIUM ST in the fram e o f ESA ITT for Robust LPV for lau ncher 
application, see lASTRIUM Sll ri2009allbh: [Ganet-Schoeller and Maurice! (l201lh. and 



for rob ust stability analysis for ATV industrial validation, see iGanet-Schoeller et al. 
(l2009h : 



results by NGC on rendez v ous va lid ation proces s in th e frame of the ESA V VAF 
project, see lDi Sotto et al.l (l2010al]bh: iKron et all (IpOloh: [Paulino et ail (boioh. and 
on re-entry validation, see IKron and de Lafontaind (120041 ) . 



Many other non industrial applications are available in the literature. Two of them 
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could be consider ed as major for our applications: a PhD thesis on miss i le app lication by 
Adounkpg (120041 ) and the IQC toolbox developed by IScherer and Kosd (120081) with Delf t 
University that was applied for HL-20 re-entry vehicle analysis in lVeenman et al.l (120091 ). 
under ESA contract. 

In this paper we propose an original approach to validation of control laws based on 
measures and convex optimization over linear matrix inequalities (LMIs). The approach 
can be seen as a blend between worst case validation techniques and statistical simulation 
techniques. On the one hand, we optimize over worst case admissible trajectories to 
validate a system property. On the other hand, we propagate along system trajectories 
statistical information (probability measures) instead of deterministic initial conditions. 
Moreover, our approach is primal in the sense that we optimize directly over systems 
trajectories, we do not seek a dual Lyapunov certificate. 

First we rephrase our validation problem as a robustness analysis problem, and then as 
a nonconvex nonlinear op timization problems over a dmissible trajectories. This is the 



approach followed e.g. in iPrajna and Rantzerl (120071 ) where the authors verify or prove 



temporal properties such as safety (all trajectories starting from a set of initial conditions 
never reach a set of bad states), avoidance (at least one trajectory starting from a set of 
initial conditions will never reach a set of bad states), eventuality (all trajectories starting 
from a set of initial conditions will reach a set of good states in finite time) and reachability 
(at least one trajectory starting from a set of initial conditions will reach a set of good 
states in finite time). 



Following lLasserre et al.l ( 120081 ) . we formulate our nonconvex nonlinear trajectory opti- 
mization problem as a linear programming (LP) problem on measures, with at most three 
unknowns: the initial measure (modeling initial conditions), the occupation measure (en- 
coding system trajectories) and the terminal measure (modeling terminal conditions); in 
some cases a subset of these decision variables may be given. The final time can be 
given or free, finite or infinite. If all the data are polynomials (performance measure, 
dynamics, constraints), then we formulate the infinite-dimensional LP on measures as a 
generalized moment problem (GMP), and we approach it via an asymptotically converg- 
ing hierarchy of finite-dimensional convex LMI relaxations. These LMI relaxations are 
then modeled with our specialized software GloptiPoly 3, solved with a general-purpose 
semidefmite programming (SDP) solver, and system properties are validated from the 
numerical solutions. 

The models we can deal with should be piecewise polynomial (or rational) in time and 
space. We consider a compact region of the state-space over which the trajectories are 
optimized. The system dynamics are defined locally on explicitly given basic semialgebraic 
sets (intersections of polynomial sublevel sets) by polynomial vector fields. The objective 
function (stability or performance measure) consists of a polynomial terminal term and 
of the time integral of a polynomial integrand. 



Our contribution is to use the approach desc ribed in lLasserre et al.l (120081 ) jointly with the 
corresponding software tool GloptiPoly 3 by Henrion et al.l ( 20091 ) to address temporal val- 
idation/verification of control systems properties, in the applied conte xt of space launcher 
applic ations. The use of LMI and measures was already investigated in lPrajna and Rantzer 
(120071 ) :'or building Lyapunov barrier certificates, and based on a dual to Lyapunov's the- 
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orem described in iRanzterl (120011 ) . Our approach is similar, in the sense that optimization 
over systems trajectories is formulated as an LP in the infinite-dimensional space of mea- 
sures. This LP problem is then approached as a generalized moment problem via a 
hierarchy of LMI relaxations. 

Historically, the idea of reformulating nonconvex nonlinear ordinary differential equations 
(ODE) into convex LP, and especially linear partial differential equations (PDE) in the 
space of probability measures, is not new. It was Joseph Liouville in 1838 who first intro- 
duced the linear PDE involving the Jacobian o f the t ransf ormation exer t ed by the so lution 



of an ODE on its initial condition iLiouville I ( 118381 ). see lEhrendorferl (119941 120021 ) for a 



survey on the Liouville PDE with applications in meteorology. The idea was then largely 
expanded in Hen ri Poinca r e's w ork on dynamical systems at the end of the 19th century, 
see in particular (|Poincare1 . ll899l. Chapitre XII (Invariants inte gr aux)). This work was pur- 
sued in the 20th century in iKrvloff and Bogolioubofil (119371 ). (INemystkii and Stepanovl . 
19471 Chapter VI (Systems with a n integral invariant)) and more re cently i n the c on- 



text of optimal trans port by e.g. iRachev an d Riis chendorJ (119981 ). IVillanil ( 120031 ) or 



Ambrosio et al. 



(2008). Poincare himself in (jPoincarel Il908l . Section IV) mentions the 
potential of formulating nonline ar ODEs a s linea r PDEs, a nd this programme h a s bee n 
carried out to some exte n t bv ICarlemanl (|1932l). see also lLasota and Mackevl (ll985|). 



Kowalski and Steebl (119911) . iHernandez-Lerma and Lasserrd (12003 ) and more recently |Baxkley_et_al 
(120061 ) . IVaidya and Mehtal (120081 ). iGaitsgory and Quincampoix! (120091). For recen t stud- 



ies of the Lio u yille equation in optimal control see e.g. iKwee and Schmidhuberl ( 120021 ) 
and iBrockettl (120071). and for application s in u ncert ainty propagation and control l aw 
validation see e.g. iMellodge and Kachrool ( 120081 ) and lHalder and Bhattacharyal (120 111 ). 



2 Piecewise polynomial dynamic optimization 

Consider the following dynamic optimization problem with piecewise polynomial differ- 
ential constraints 

J = inf x (t) h T (T,x(T)) + J Q T h(t,x(t))dt 

s.t. x(t) = f k (t,x{t)), x{t) eX k , k = 1,2,..., N (1) 

x(0) g x , x(T) el T , te [o,T] 

with given polynomial dynamics f k G x] and costs h, hx G M.[t, x], and state trajectory 
x(t) constrained in closed basic semialgebraic sets 

X k = {x e R n : g kj (t,x) > 0, j = 1,2,.. .,N k } 

for given polynomials g k j G M.\t, x]. We assume that the state-space partitioning sets, or 
cells X k , are disjoint, i.e. all their respective intersections have zero Lebesgue measure in 
R™, and they all belong to a given compact semialgebraic set, e.g. 

X = {xeR n : \\x\\l < M} 

for a sufficiently large constant M > 0. Finally, initial and terminal states are constrained 
in semialgebraic ets 

X = {x G R n : g 0j (t, x) > 0, j = 1, 2, . . . , N } c X 
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and 

X T = {xeR n : g Tj (t,x)>0, j = 1,2,...,N T } CX 
for given polynomials g j,gTj G M[t,x]. 

It is assumed that optimization problem ([T]) arises from validation of a controlled (closed- 
loop) system behavior. In this problem, the infimum is sought over absolutely continuous 
state trajectories x(t), but the same methodology can be extended to (possibly discon- 
tinuous) trajecto ries of bound e d var iation (e.g. in the presence of impulsive controls or 
state jumps), see lClaeys et al.l (120111 ). 

The final time T is either given, or free, in which case it becomes a decision variable, 
jointly with x(t). Similarly, the initial and terminal constraint sets X and X? are either 
given, or free, in which case they also become optimized decision variables. 



3 Linear programming on measures 



In this section, we formulate optimization problem ([[]), which is nonlinear and nonconvex 
in state trajectory x, into a convex optimization problem linear in /i, the occupation 
measure of trajectory x. This reformulation is classical in the m odern theory o f calculus 
of variations and op ti mal control a nd it can be traced back to lYoungi (119371 ). see also 
Ghouila-Houril (119671 ). lYoungi (119691 ) and iGamkrelidzd (119751 ). amongst many others. 



3.1 Occupation measure 

To understand the basic idea behind the transformation and the elementary concept of 
occupation measure, it is better to deal with the single nonlinear ODE 



it = f(xt) 



(2) 



where Xt is a shortcut for x(t), a time-dependent state vector of R n and / : R™ — >■ R n is a 
uniformly Lipschitz map. It follows that the Cauchy problem for ODE (jSJ) has a unique 
solution x t for any given initial condition i € K". 

Now think of initial condition Xq clS Sb random variable of R n , or more abstractly as a 
nonnegative probability measure /xq with support X C R n , that is a map from the sigma- 
algebra of subsets of X to the interval [0, 1] C R such that [Iq(Xq) = 1. For example, the 
expected value of xo is the vector E[xq] = j x xno{dx). 

Now solve ODE (j2J) for a trajectory, or flow x t , given this random initial condition. At 
each time t, state Xt can also be interpreted as a random variable, i.e. a probability 
measure that we denote by ji t {dx). We say that the measure is transported by the flow 
of the ODE. The one- dimensional family, or path of measures fit satisfies a PDE 



^ + div(M) = o 



(3) 



which turns out to be linear in the space of pro bability meas ures. This PDE is usually 
called Liouville's equation. As explained e.g. in (jVillanil . I2003L Theorem 5.34), nonlinear 
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ODE (JT]) follo ws by applyin g Cauchy's method of characteristics to linear transport PDE 



see also (jEvansl . l2010l Section 3.2) for a tutorial exposition. In equation (j3J), div 



or distributional sense, i.e. 



models the divergence operator, i.e. 

n g 

div W = E £: 

for every smooth function v. Its action on measures should be understood in the weak, 

J v div(z/) = — J Dv ■ dv 

where v is a smooth test function, v is a vector-valued measure and D stands for the 
gradient operator. Given a subset T X X in the sigma- algebra of subsets of [0, T] x X, 
we define the occupation measure 

|i(Tx Af) = J [i t {X)dt 

which encodes the time-space trajectories x t , in the sense that /x([0,T] x X) is the total 
time spent by trajectory x t in a subset X C X of the state space. In its integral form, 
transport PDE (J3]) becomes 

div(/» = /x - Mr (4) 

where /x^ is the terminal probability measure with support Xr C R n . PDE (J3J) can 
equivalently be formulated as 



Dt; • /d/x = / vdfiT — / t>d/x (5) 

for all smooth functions v compactly supported on X. Problem ([5]) is an infinite-dimensional 
linear system of equations linking occupation measure /x, initial measure /x and terminal 
measure /xt, consistently with ODE (J2]). 



3.2 Finite terminal time 

If we apply these ideas to problem (TJp), encoding the state trajectory a;(i) in an occupation 
measure /x supported on X, we come up with an infinite-dimensional LP problem 

Joo = inf M / h T (T, x)dfi T (x) + Efc J ' x )d^k(t, x) 

/ £>f ar) ■ /*(*, a;) d/x*.(t, a;) = 
J v(T,x) d[Pr(x) — J v(0, x) dfj, (x) 

for all smooth test functions v, where each occupation measure /Xfc is supported on set X^ 
and the global occupation measure is 

A 4 — E ^ k 

k 
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with normalization constraint 

M[0,T]xX)=T (7) 
such that T is a finite terminal time. 

In problem ([6]), final time T, initial measure /io and terminal measure \xt may be given, 
or unknown, depending on the original optimization problem. 



3.3 Generalized moment problem 

More concisely, LP problem ([6]) can be formulated as follows: 

S/t - Sfc Ix k aki d ^ k = ^ V ^ 

where the unknowns are a finite set of nonnegative measures with respective compact 
semialgebraic supports 

X k = {x : g kj (x) > 0, Vj}. (9) 

Note that here we have incorporated time variable t into vector x, for notational concise- 
ness. If all the data in problem (|6]) are polynomials, and if we generate test functions 
v{x) using a polynomial basis (e.g. monomials, which are dense in the set of continuous 
functions with compact support), all the coefficients a(x), b(x), c(x) are polynomials, and 
there is an infinite but countable number of linear constraints indexed by i. 

We will then manipulate each measure //& via its moments 

Vka = / x a d/i h (x), Va (10) 



x 



gathered into an infinite-dimensional sequence y k indexed by a vector of integers a, where 
we use the multi-index notation LP measure problem ([H]) becomes an LP 

moment problem, or GMP, see lLasserrel ( 20091 ): 



provided we can handle the repr e senta tion condition (flOj) which links a measure with its 



moments. It turns out (ILasserrd . 120091 . Chapter 3) that if sets X k are compact semialge- 
braic as in ([9]), we can use results of functional analysis and real algebraic geometry to 
design a hierarchy of LMI relaxations which is asymptotically equivalent to the generalized 
moment problem. In each LMI relaxation 

J d = infy c T y 

s.t. Ay = b 

M d (y)h0 

M d (g kj ,y)h0, Vj,k 

we truncate the moment sequence to a finite number of its moments. If the highest 
moment power is 2d, we call LMI relaxation of order d the resulting finite-dimensional 
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truncation. Matrix M d (y) is symmetric and linear in y, it is called the moment matrix, 
and it must be positive semi definite for (fTUj) to hold. Symmetric matrices M d (g k j,y) are 
also linear in y and they are called localizing matrices. They ensure that the moments 
correspond to measures with appropriate supports. For d finite, problem (1TT]1 is a finite- 
dimensional convex LP problem in the cone of positive semidefinite matrices, or SDP 
problem, for which off-the-shelf solvers are avail able, in particular implementation s of the 
primal-dual interior-point methods described in iNesterov and Nemirovskiil (11994J ). 



3.4 Convergence 

The hierarchy of LMI relaxations ( I lip generates an asymptotically converging monoton- 
ically increasing sequence of lower bounds on LP flS]), i.e. J d < J d+1 for all d = 1, 2 . . . 
and lim^oo J d = Joo- Obviously, the relaxed LP cost in problem ^ is a lower bound on 
the original cost in problem (pQ), i.e. Jqo < J. Under mild assumptions, we can show that 
actually = J, but the theoretical background required to prove this lies beyond the 
scope of this application paper. 



3.5 Long range behavior 



In the case terminal time T tends to infinity, normalization constraint (j7|) becomes irrele- 
vant, and the overall mass of occupation measure \x tends to infinity. A more appropriate 
formulation consists then of normalizing measure \x to a probability measure 



7T 



T 



k 



so that PDE (HI) becomes 



div(/vr) 



T^oo T 



and LP problem (jH]) becomes 



s-t. J2k I ^af 1 d7l k(x) 



Efe J -^dTr k (x) + 
J2 k j Dv(t, x) ■ f k (t, x) dn k (t, x) 



where test functions v(t, x) satisfy appropriate integrability and/or periodicity properties. 
Measure 7r is called an invariant probability measure and it encodes equilibrium points, 

period i c orbits, ergodic behavio r , poss i bly chaotic attractors, see e.g. |Lasota and Mackev 

(11985 ). Diaconis and Freedman ( 1999 ). Hernandez- Lerma and Lasserrel ( 2003 ) and Gaitsgorv and Quin 
(l2009h . 
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4 Application to orbital launcher control law valida- 
tion 



In this section we report the application of our method to simplified models of an atti- 
tude control system (ACS) for a launcher in exo-atmospheric phase. We consider first a 
one degree-o f-freedom (1D0 F) model, and then a 3D0F model. The full benchmark is 
described in (ISAFEVl . l201lL Section 9) but it is not publicly available, and our computer 
codes cannot be distributed either. 



4.1 ACS 1DOF 

First we consider the simplified 1DOF model of launcher ACS in orbital phase. The closed- 
loop system must follow a given piecewise linear angular velocity profile. The original 
benchmark includes a time-delay and a pulsation width modulator in the control loop, 
but in our incremental approach we just discard these elements in a first approximation. 
See below for a description of possible extensions of our approach to time-delay systems. 
As to the pulsation width modulator, it can be handled if appropriately modeled as a 
parametric uncertainty. 

4.1.1 Model 

The system is modeled double integrator 

I6(t) = u{t) 

where J is a given constant inertia and u(t) is the torque control. We denote 

e(t) 



x(t) 

;le xi(i 

the torque control is given by 



6(t) 

and we assume that both angle x\(t) and angular velocity x 2 {t) are measured, and that 



u(x(t)) = sat(K T dz(x r (t) - x(t))) 

where x r (t) is the reference signal, K e IR 2 is a given state feedback, sat is a saturation 
function such that sat(y) = y if \y\ < L and sat(y) = L sign(y) otherwise, dz is a dead- 
zone function such that dz(x) — if \xi\ < Di for some 2 = 1,2 and dz(x) = 1 otherwise. 
Thresholds L > 0, D\ > and D 2 > are given. 

We would like to verify whether the system state x(t) reaches a given subset Xt = {x G 
R 2 : x T x < e} of the deadzone region after a fixed time T, and for all possible initial 
conditions x(0) chosen in a given subset X of the state-space, and for zero reference 
signals. We formulate an optimization problem ([T]) with systems dynamics defined as 
locally affine functions in three cells Xk, k = 1,2,3 corresponding respectively to the 
linear regime of the torque saturation 



Ii = {x6R 2 : \K 2 x\ < L}, 



Xi 

-K T x 



9 



the upper saturation regime 



X 2 = {x E K 2 : K T x > L}, f 2 (x) = 

and the lower saturation regime 

X 3 = {x E K 2 : K T x < -L}, f 3 (x) = 

The objective function has no integral term and a concave quadratic terminal term 
hx{x) = —x(T) T x(T) which we would like to minimize, so as to find trajectories with 
terminal states of largest norm. If we can certify that for every initial state x(0) chosen in 
X the final state x(T) belongs to set included in the deadzone region, we have validated 
our controlled system. 




L 



4.1.2 Validation script 



The resulting GloptiPoly 3 script, implementing some elementary scaling strategies to 
improve numerical behavior of the SDP solver, is as follows: 



I = 27500; °/„ inertia 

kp = 2475; kd = 19800; % controller gains 
L = 380; % input saturation level 

dzl = 0.2*pi/180; dz2 = 0.05*pi/180; % deadzone levels 
thetamax = 50; omegamax =5; % bounds on initial conditions 
epsilon = sqrt(le-5); % bound on norm of terminal condition 
T = 50; % final time 



d = input ('order of relaxation ='); d = 2*d; 



% measures 
mpoK'xl' ,2) 
mpol('x2' ,2) 
mpol('x3' ,2) 
mpoK'xO' ,2) 
mpoK'xT' ,2) 



ml = meas(xl) 
m2 = meas(x2) 
m3 = meas(x3) 
mO = meas(xO) 
mT = meas(xT) 



7o linear regime 
7o upper sat 
% lower sat 
7o initial 
7o terminal 



7o dynamics on normalized time range [0,1] 
7o saturation input y normalized in [-1,1] 
K = -[kp kd]/L; 

yl = K*xl; fl = T*[xl(2); L*yl/I] ; 7. linear regime 
y2 = K*x2; f2 = T*[x2(2); L/I] ; % upper sat 
y3 = K*x3; f3 = T*[x3(2); -L/I]; % lower set 

7o test functions for each measure = monomials 

gl = mmon(xl,d); g2 = mmon(x2,d); g3 = mmon(x3,d); 

gO = mmon(x0,d); gT = mmon(xT,d) ; 
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7, unknown moments of initial measure 
yO = mom(gO) ; 

°/ unknown moments of terminal measure 
yT = mom(gT) ; 

% input LMI moment problem 

cost = mom(xT'*xT); 

Ay = mom(diff (gl,xl)*f 1)+. . . 

mom(diff (g2,x2)*f2)+. . . 

mom(diff (g3,x3)*f3) ; % dynamics 
°/ trajectory constraints 
X = [yl"2<=l; y2>=l; y3<=-l] ; 
°/o initial constraints 

XO = [xO(l) ~2<=thetamax~2, x0(2) ~2<=omegamax~2] ; 

% terminal constraints 

XT = [xT , *xT<=epsilon~2] ; 

°/o bounds on trajectory 

B = [xl'*xl<=4; x2'*x2<=4; x3 J *x3<=4]; 

°/ input LMI moment problem 
P = msdp (max (cost) , ... 

mass(ml)+mass(m2)+mass(m3)==l, . . . 

mass(mO)==l , ... 

Ay==yT-yO, ... 

X, XO, XT, B); 

% solve LMI moment problem 
[status, obj] = msol(P) 



4.1.3 Numerical results 

With this Matlab code and the SDP solver SeDuMi 1.3 (see sedumi.ie.lehigh.edu) 
we obtain the following sequence of upper bounds (since we maximize) on the maximum 
squared Euclidean norm of the final state: 



relaxation order d 


1 


2 


3 


4 


upper bound Jd 


1.0 • 10~ 5 


1.0 • lO" 5 


1.0 • HT 5 


1.0 ■ 10~ 5 


CPU time (sec.) 


0.2 


0.5 


0.7 


0.9 


number of moments 


30 


75 


140 


225 



In the table we also indicate the CPU time (in seconds, on a standard desktop computer) 
and the total number of moments (size of vector y in the LMI relaxation ( II ip . We see that 
the bound obtained at the first relaxation (d = 1) is not modified for higher relaxations. 
This clearly indicates that all initial conditions are captured in the deadzone region at 
time T, which is the box [-2,2]±^ x [-5,5]^^ D {i 6 R 2 : x T x < 10~ 5 }. 



11 



If we want to use this approach to simulate a particular trajectory, in the code we must 
modify the definition of the initial measure. For example for initial conditions £i(0) = 50, 
£2(0) = —1, we must insert the following sequence: 

% given moments of initial measure = Dirac at xO 
p = genpow(3,d); p = p(:,2:end); % powers 
thetaO = 50; omegaO = -1; % in degrees 
yO = ones (size (p , 1) , 1)* [thetaO omegaO] *pi/180 ; 
yO = prod(y0.~p,2) ; 

As previously, the sequence of bounds on the maximum squared Euclidean norm of the 
final state is constantly equal to 1.0 • 10 -5 , and in the following table we represent as 
functions of the relaxation order d the masses of measures fi^, k = 1,2,3 which are 
indicators of the time spent by the trajectory in the respective linear, upper saturation 
and lower saturation regimes: 



relaxation order d 


1 


2 


3 


4 


5 


6 


7 


J dfj,i 


37 


89 


92 


92 


93 


93 


93 


Jdfi 2 


32 


5.3 


0.74 


0.30 


0.21 


0.15 


0.17 


JdtM 3 


32 


5.1 


7.1 


6.9 


6.8 


6.9 


7.0 



This indicates that most of the time (approx. 93%) is spent in the linear regime, with 
approx. 7% of the time spent in the lower saturation regime, and a negligible amount of 
time is spent in the upper saturation regime. This is confirmed by simulation, see Figure 

m 

4.2 Dealing with uncertainty 

Finally, note that we can incorporate real parametric uncertainty in the dynamics, if 
required. Each uncertain parameter must be introduced as an additional state of the 
system, and this generally makes the dynamics polynomial in the extended state. The 
parameters must be constrained to an explicitly given compact semialgebraic set, and 
the overall trajectory optimization problem consists of finding the worst-case uncertain 
parameter instance. 

For illustration, we modify the previous script to cope with uncertainty entering affmely 
the dynamics. We assume that the (reciprocal of the) inertia is subject to multiplicative 
relative uncertainty, i.e. we replace occurences of / with jt^I for u a real parameter such 
that u 2 < U with U > a given threshold, say U = |. The resulting GloptiPoly 3 script 
is as follows: 



I = 27500; % inertia 

kp = 2475; kd = 19800; % controller gains 
L = 380; % input saturation level 

dzl = 0.2*pi/180; dz2 = 0.05*pi/180; % deadzone levels 
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_400 1 1 1 1 1 1 1 1 1 1 1 

5 10 15 20 25 30 35 40 45 50 

time 

Figure 1: Torque input with lower saturation during approx. 7% of the time range. 
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thetamax = 50; omegamax =5; % bounds on initial conditions 
epsilon = sqrt(le-5); % bound on norm of terminal condition 
T = 50; % final time 

ui = 0.5; % uncertainty level (in percent) 

d = input('order of relaxation ='); d = 2*d; 

°/ measures = states + uncertain parameter 
mpol( 'xl ' ,2) ; mpol ul; ml = meas (xl ,ul) ; °/„ linear regime 
mpol( 'x2' ,2) ; mpol u2; m2 = meas(x2,u2); % upper sat 
mpol( 'x3' ,2) ; mpol u3; m3 = meas(x3,u3); % lower sat 
mpol( 'x0' ,2) ; mpol uO; mO = meas(x0,u0); °/„ initial 
mpol( 'xT' ,2) ; mpol uT; mT = meas(xT,uT); °/ terminal 

% dynamics on normalized time range [0,1] 
°/ saturation input y normalized in [-1,1] 
K = -[kp kd]/L; 

yl = K*xl; fl = T*[xl(2); (1+ul) *L*yl/I] ; °/„ linear regime 
y2 = K*x2; f2 = T*[x2(2); (l+u2)*L/I] ; % upper sat 
y3 = K*x3; f3 = T*[x3(2); -(l+u3)*L/I] ; % lower set 

°/ test functions for each measure = monomials 

gl = mmon( [xl ;ul] ,d) ; g2 = mmon( [x2;u2] ,d) ; g3 = mmon( [x3;u3] ,d) ; 
gO = mmon( [x0;u0] ,d) ; gT = mmon( [xT;uT] ,d) ; 

°/ unknown moments of initial measure 
yO = mom(gO) ; 

°/ unknown moments of terminal measure 
yT = mom(gT); 

°/ input LMI moment problem 

cost = mom(xT'*xT); 

Ay = mom(diff (gl,xl)*f 1)+. . . 

mom(diff (g2,x2)*f2)+. . . 

mom(diff (g3,x3)*f3) ; % dynamics 
°/ trajectory constraints 
X = [yl-2<=l; y2>=l; y3<=-l] ; 
°/ initial constraints 

X0 = [x0(l) "2<=thetamax"2, x0(2) "2<=omegamax"2] ; 

% terminal constraints 

XT = [xT , *xT<=epsilon~2] ; 

°/ bounds on trajectory 

B = [xl'*xl<=4; x2'*x2<=4; x3 J *x3<=4]; 

% bounds on uncertain parameter 

U = [ul~2<=ui~2; u2~2<=ui~2; u3"2<=ui"2] ; 

UO = [u0"2<=ui"2] ; UT = [uT~2<=ui~2] ; 
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°/ input LMI moment problem 
P = msdp (max (cost) , ... 

mass(ml)+mass(m2)+mass(m3)==l, . . . 

mass(mO)==l , ... 

Ay==yT-yO, ... 

X, XO, XT, B, ... 
U, UO, UT) ; 

% solve LMI moment problem 
[status, obj] = msol(P) 



Introducing an uncertain parameter amounts to adding a state variable to the model, and 
this has a significant impact on the overall computational time, as shown in the table 
below: 



relaxation order d 


1 


2 


3 


4 


upper bound Jd 


1.0- 1(T 5 


1.0 • 10" 5 


1.0 • KT 5 


1.0 ■ 10~ 5 


CPU time (sec.) 


0.8 


2.4 


7.3 


17.7 


number of moments 


75 


224 


501 


946 



Note however that for this example there is no effect on the bounds, and the control law 
is also validated in the presence of uncertainty on inertia. 



4.3 ACS 3DOF 

The main difficulty with the 3 degree-of-freedom version of the ACS benchmark is the 
large number of states (4 quaternions for the launcher attitude, 3 angular velocities, and 
2 additional states for the quaternion reference signal) and the non-linear (quadratic) 
dynamics. 



4.3.1 Model 

The following example illustrates the validation of a given control law following a given 
constant roll velocity equal to u>f = 20°/s, in the absence of actuator saturation, but 
in the presence of non-linear coupling between the three axes of the launcher. System 
dynamics are given by x = f(x) with 





1 




fq( X ) 


X = 


w 


, f(x) = 


fw(x) 




q R 




_ fq*(x) _ 



and 

f q (x) = ^Q(x) w , 
f w (x) = I^iutyx) - Q(x)Iw), 
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where 



Q(x) 



f q 4 x ) = 
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-92 
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92 
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93 
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9i 


9o 
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gf 
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and g e ffi 4 is the quaternion modeling the attitude, wsK 3 models the angular velocities, 
is the reference quaternion to follow, which depends on the constant reference 
is the quaternion error, and w E G M 3 is the angular velocity 



q R e 

velocity u>f G 

error. Both vectors q(t) and q R (t) have constant Euclidean norm since J7II9II2 — 2q T q = 



q E G 



q T f q (x) = for all x and similarly 



A 

dt 



n R\\ 2 
9 II2 



0, so we enforce the algebraic constraints 



hi; 
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The proportional-derivative control law u(x) is designed so that velocity Wi(t) follows 
reference velocity , and our validation tasks consists of optimizing over the worst-case 
trajectory starting at a given initial condition x(0) and maximizing the objective function 



J 



T 



(wi(t) - w^fdt 



for a given time horizon T. If we can guarantee a sufficiently small upper bound Jd on this 
objective function, for a given relaxation order d — 1, 2, . . ., the control law is validated. 



4.3.2 Validation script 

We use the following GloptiPoly 3 script (note however that the function G2I to generate 
the inertia data is not available for the reader): 

% Inertia 

Ig = [ 27500, -50, -1100; 

-50, 45300, -220; 
-1100, -220, 44100] ; 
[Ge2In,Ip] = G2I(Ig); 

% controller gains 
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wO = 0.3; xi = 0.7; 

Kpl = Ip(l,l)*w0~2; Kdl = 2*Ip(l , 1) *xi*wO; 
wO = 0.25; xi = 3; 

Kp2 = Ip(2,2)*w0~2; Kd2 = 2*Ip (2 , 2) *xi*wO ; 

Kp3 = Ip(3,3)*w0~2; Kd3 = 2*Ip (3 , 3) *xi*wO ; 

°/o initial conditions 

QO = [1;0;0;0]; % attitude quaternion 

WO = [0;0;0] .* (pi/180) ; 7. angular velocities 

T = 50; % final time 

% reference velocity 
WR1 = 20*pi/180; 

d = input ('order of relaxation ='); 
°/ occupation measure 

mpol('q',4); mpoK'w',3); mpoK'qr' ,2) ; 
x = [q;w;qr]; m = meas(x); 

°/ initial measure 

mpoK'qO' ,4) ; mpoK'wO' ,3) ; mpol( 'qrO' ,2) ; 
xO = [qO;wO;qrO]; mO = meas(xO); 

°/ terminal measure 

mpoK'qt' ,4) ; mpoK'wt' ,3) ; mpol( 'qrt ' , 2) ; 
xt = [qt;wt;qrt]; mt = meas(xt); 

°/ dynamics on normalized time range [0,1] 
Omega = [0 -w(3) w(2) ; w(3) -w(l) ; -w(2) w(l) 
Qr = [qr(l) qr(2) 0; 

-qr(2) qr(l) 0; 

qr(l) qr(2); 

-qr(2) qr(l)] ; 
dq = 2*Qr*q; 

Q = [q(l) -q(2) -q(3) -q(4) ; 

q(2) q(l) -q(4) q(3) ; 

q(3) q(4) q(l) -q(2) ; 

q(4) -q(3) q(2) q(l)] ; 
fq = 1/2*Q* [0;w] ; 
Torque = [-Kdl* (w(l) -WR1) ; 

-Kd2*(w(2)-0)-Kp2*2*dq(3) ; 
-Kd3*(w(3)-0)-Kp3*2*dq(4)] ; 
fw = Ig\ (Torque-Omega*Ig*w) ; 
fqr = [-l/2*qr(2)*WRl;l/2*qr(l)*WRl] ; 
f = T*[fq;fw;fqr] ; 

°/ test functions 
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g = mmon(x,d) ; 
gO = mmon(xO,d) ; 
gt = mmon(xt ,d) ; 

% given moments of initial measure (Dirac) 
p = genpow(10,d) ; p = p(:,2:end); % powers 
yO = ones(size(p,l) ,1)*[Q0 ) WO' 1 0]; 
yO = prod(y0.~p,2); 

°/ unknown moments of final measure 
p = genpow(10,d) ; p = p(:,2:end); % powers 
yt = ones (size (p, 1) , 1) * [qt ' wt' qrt']; 
yt = mom (prod ( (yt . ~p) ' ) ' ) ; 

% input LMI moment problem 

cost = mom( (w(l) -WR1) ~2) ; % objective function 

Ay = mom(dif f (g,x) *f ) ; % dynamics 

X = [q'*q == 1; qr'*qr == 1]; % quaternion 

Xt = [qt J *qt == 1; qrt ' *qrt == 1]; % final quaternion 

P = msdp (max (cost) , mass(m)==l, yt==Ay+yO, X, Xt) ; 

°/ solve LMI moment problem 
[status, obj] = msol(P) 

4.3.3 Numerical results 

We obtain the monotonically decreasing upper bounds Jd, d — 1, . . . ,4 reported in the 
following table: 



The first LMI relaxation (110 moments of degree up to 2) seems to be unbounded above 
(the dual LMI problem is infeasible) and hence it conveys no useful information. The 
second LMI relaxation (770 moments of degree up to 3) and the third LMI relaxation 
(1430 moments of degree up to 4) are solved with SeDuMi 1.3 in a few seconds. The 
fourth LMI relaxation is more challenging, with 5720 moments of degree up to 5, and it 
requires less than one hour of CPU time on a standard PC. We observe however that a 
useful upper bound on J is already obtained at a low relaxation order (say d = 2 or d = 3) 
at a low computational cost, and that the computational burden increases significantly 
for d — 4 without dramatic improvements in the quality of the upper bound. 



relaxation order d 



12 3 4 



upper bound Jd 
CPU time (sec.) 
number of moments 



oo 7.5198 -10- 3 2.9110 -10- 3 2.9085 • 10" 3 
0.2 11.6 24.2 2640 

110 770 1430 5720 
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4.3.4 Coping with saturations and dead-zones 

We can modify the above code to cope with actuator saturation and dead- zone, proceeding 
exactly as described above for the ACS 1D0F benchmark problem, namely by splitting 
the trajectory occupation measures into local occupation measures defined on semialge- 
braic (here polyhedral) cells. In the presence of saturations on the three actuators, this 
generates 3 3 = 27 local occupation measures. This should not be an issue from the com- 
putational point of view, the limiting factor being essentially the size of the largest SDP 
block, which grows polynomially with the relaxation order, but with an exponent which 
is the number of variables (here equal to 10). 

4.3.5 Following a time-varying reference signal 

Finally, if we want to validate a control law to follow a time-varying velocity profile 
(instead of a constant velocity as above), we must introduce time as an additional variable 
in the trajectory occupation measure, we must introduce the velocity as an additional 
state, and we must split the occupation measure into local occupation measures, each of 
which corresponding to given time-invariant dynamics. Alternatively, we can also model 
time-varying dynamics, but the dependence on time must be polynomial (which is not 
the case e.g. if the reference signal is piecewise linear). 

5 Discussion 

In this section we discuss the limitations of our approach and we also describe possible 
extensions. 

5.1 Limitations 

The main limiting factor for our method is the total number of variables (number of 
states plus number of uncertain parameters), since the computational burden of solving 
the LMI relaxations grows polynomially as a function of the LMI relaxation order, but 
the order of the polynomial dependence depends linearly on the number of variables in 
the measures (here the number of states). So the critical dimension is not too much 
the number of measures (that is the number of cells partioning the state-space, which 
corresponds to our models of nonlinearities), or the degree of the polynomials in the cost 
and/or dynamics, but rather the number of states. Systems with a high number of states 
can be handled with these techniques, but only when exploiting problem structure and 
sparsity. 

More explicitly, a rough complexity analysis can be carried out as follows. If an LMI 
relaxation has the simplified form 

inf y c T y 

s.t. M d (y) y 
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where M<i(y) is the moment matrix of a measure of n variables at relaxation order d, then 
the number of variables in vector y is 



N = 

and the size of matrix Md(y) is 

M 



n + 2d\ _ (n + 2d)\ 
n ) ~ n \ (2d)! 




(n 



+ d)\ 



n 



\d\ 



A standard primal-dual interior-point algorithm to solve this LMI at given relative accu- 
racy e > (duality gap threshold) requires a number of iterations (Newton steps) growing 
as 0(M2 loge), whose dependence on M is sublinear, hence almost negligible. In prac- 
tice, on well-condition ed problems, we observe that th e number of Newton iterat i ons is 



bet ween 5 and 50. seelNesterov and Nemirovskiil ( 119941 ) , IVandenberghe and Boydl (119961 ) 
and Ben-Tal and Nemirovski ( 200ll ). 



In the real model of computation (for which each addition, subtraction, multiplication, di- 
vision of real numbers has unit cost), each Newton iteration requires 0(N 2 M)+0(N 3 M) + 
0(N 2 M 2 ) operations to form the linear system of equations, and 0(M 3 ) operations to 
solve the system and find the search direction. When solving a hierarchy of simple LMI 
relaxations as described above, the number of variables n is fixed, and the relaxation order 
d varies, so the dominating term in the complexity estimate grows in 0(d 4n ), which clearly 
shows a strong dependence on the number of variables. Even though the growth of the 
computational burden is polynomial in the relaxation order, the exponent is 4 times the 
number of variables. One should however keep in mind that these estimates are (usually 
very loose) asymptotic upper bounds, and that the observed computational complexity 
grows much more moderately in practice (at least in the absence of conditioning and 
numerical stability issues). 

Finally, if the LMI relaxation has the form 

iniy c T y 

s.t. M d {y k )h0, k = l,2,...,K 

where M^yu) is the moment matrix of a measure fi^ of n variables at relaxation d, and 
we have K measures, then the above complexity estimate grows in 0(Kd An ). The impact 
on the computation burden of the number n of variables in each measure fik is thus much 
more critical than the number K of measures. In our target application, n is the number 
of states and uncertain parameters, K is the number of cells used to model nonlinearities, 
and a lower bound on the minimum relaxation order d is given by the degree of the 
polynomial data (dynamics, constraints, objective function). 

5.1.1 Accuracy 

The original validation problem can be formulated as infinite-dimensional linear program- 
ming problems on the dual cones of nonnegative measures and continuous functions, so 
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that there is no hope of finite convergence of finite-dimensional LMI optimization tech- 
niques. Convergence is guaranteed only asymptotically. Obviously, the accuracy depends 
on the speed of convergence of the hierarchy of LMI problems, but this speed is impossible 
to evaluate a priori. The only guarantee is that we have a monotically increasing sequence 
of lower bounds on the objective function to be minimized. 

As far as solving finite-dimensional LMI problems is concerned, we should emphasize the 
fact that there is currently no proof of backward stability of implementations of LMI 
solvers. Moreover, evaluation of the conditioning of a given LMI problem is difficult, and 
estimates can be obtained only at the price of solving several instances of (slightly modified 
versions) the original LMI problem. To certify the output of a numerical algorithm, we 
must simultaneously ensure backward stability of the algorithm and well-conditioning of 
the problem. None of these two properties can be ensured when solving LMI problems, 
given the current state-of-the-art in LMI solvers. But this is not specific to our method, 
and any validation-verification technique which is based on LMI optimization is subject 
to the same limitations. 



5.1.2 Time-delay systems 

Our techniques can in principle cope with time-delays, but this requires a significant 
modification of the approach described in this document. Let us only sketch the main 
ideas, in the case of an ordinary differential equation with one time-delay 

x(t) = f(x(t))+g(x(t-r)), WG[0,T] 

with boundary conditions 

x{t) = Z(t), We[-r,o] 

where £(t) is a given function recording the state history due to the given delay rEl. 
Instead of transporting a probability measure fit{dx) supported on X C W 1 from initial 
time t — to terminal time t = T, we must transport the state history in an occupation 
measure /J, t (ds,dx) supported on [— r, 0] x X for t G [0, T]. 



5.1.3 Discrete-time systems 

In this paper we deal only with continuous-time systems. Discrete-time systems of the 
form 

%k+i = f(x k ) (12) 

must be handled differently. Denoting by Hk{x) the probability measure transported along 
dynamical system (lT2"j) . the discrete-time analogue of Liouville's transport equation 03) 
reads 

Vk+i(X) = / Vk(dx) = / I x (f(x))fi k (dx) (13) 

Jf-l(X) J 

where X is any subset of the sigma-algebra of state set X, and lx is the indicator function 
equal to one in X and zero outside. It follows from ( f!3|) that the moments of measure 
Hk+i can be expressed linearly as functions of moments of measure fik- Besides this 
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analogy, the resulting discrete-time generalized moment problem differs significantly from 
its continuous-time counterpart. For this reason, handling discrete-time systems would 
require an important modification of the approach described in this document. 



6 Conclusion 



In this paper we do not follow the mainst ream Lyapunov ap proach to dynamical systems 
stability and/or performance validation. iLvapunov I (118921 ) originally introduced his ap- 
proach to conclude about systems behaviour while avoiding explicit computa tions of sys- 
tem tr a jectories. In the 1990s it was com bined with LMI techniques, see e.g. iBoyd et al. 



(J1994]); iNesterov and Nemirovskiil ( 119941 ) to provide numerically quadratic certificates of 
stability and/or performance for linear and nonlinear systems. Later on, it was extended 
to more general, but still nu merical polynomial sum-of-squares (SOS) certificates, see e.g. 



Henrion and Garullil (120051 ) . By contrast, our approach is genuinely primal, in the sense 



that we directly optimize numerically over systems trajectories, using measures, moments 
and LMI techniques. Indeed, since the computation of SOS certificates is numerical any- 
way, we believe that in the context of validation, LMI techniques should rather be used to 
optimize systems trajectories, instead of optimizing indirect certificates for these trajecto- 
ries. Since we are using primal-dual SDP solvers to solve the hierarchy of LMI problems, 
solving the dual to the problem of optimizing trajectories provides anyway certificates of 
feasibility or infeasibility, i.e. of stability, instability and/or performance. 

As briefly described in §0 the measure/LMI approach can be extended easily to systems 
with real uncertainties, as soon as the uncertainties enter polynomially in the dynamics, 
and they are bo unded in explicitly given semialgebraic sets (e.g. balls or boxes). See e.g. 
(jLasserrd . I2009L Section 13.1) for more information on robust optimization. As usual, the 
number of real uncertain parameters should be kept small enough to ensure computational 
tractability. 

Note that in this paper the set of initial conditions is assumed to be given, and when the 
initial measure is unknown, its support is constrained to be included in the set of initial 
conditions. In the context of validation/ verification it could be relevant to maximize the 
size of the set of initial conditions, and we are currently investigating this problem from 
the point of view of occupation measures. 

The measure/LMI approach can readily deal with systems with piecewise polynomial 
(or rational) models, e.g. systems with input/output saturations and dead-zones. It 
can be seen as a (primal) extension of early attempts of us e of LMI /Lyapunov tech- 



niques in the context of systems with satur ation s, see e.g. iHenrion and Tarbouriech 



(119991 ) and more recently iGarulli et al.l (1201 if ) and iTarbouriech et al.l (120111 ). As briefly 
discussed in §5\ our approach can be extended without major difficulty to time-delay, 
stochastic, discrete-time systems or hybrid system dyna mics where transitions betwee n 
models are ruled e.g. by probability me asures, see e.g. iDiaconis and Freedman ( 1999 ). 



Hernandez- Lerma and Lasserrd (120031 ) or iBarklev et al.l ( 120061 ) . 
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